clear all
set more off, perm
local mydir "`1'"
global Tables 	`mydir'/BOOT/Tables
global Work   	`mydir'/BOOT/DataWork
global Raw    	`mydir'/BOOT/DataRaw

* uses same data as TS:
global Data     `mydir'/TS/DataForTables 
global WRDS  	`mydir'/WRDS
/* ********************************************************************* */
set varabbrev off
/* ********************************************************************* */
* Get data on return of the S&P between t-11 and t:
import sas $WRDS/CRSP/msp500.sas7bdat, clear case(lower)
gen datem=mofd(caldt)
format datem %tm

tsset datem
gen lnret=ln(1+sprtrn)

* for sanity check (make sure that returns in CRSP line up with dependent variable in Tables)
gen r1=0
forvalues x=0(1)11 {
	replace r1=r1+L`x'.lnret
}	
keep datem lnret r1  
tempfile tmp
save "`tmp'", replace
/* ********************************************************************* */
use datem statpers STGE2 D P using $Data/TimeSeries01, clear
merge 1:1 datem using "`tmp'"
keep if _merge==3
drop _merge
keep if dofm(datem)<=mdy(01,1,2021) 
keep if dofm(datem)>=mdy(12,1,1981)

* find "rho" to compute cumulative log returns 
gen dp=ln(D)-ln(P) if STGE2<.
quietly: su  dp 
scalar rho=1/(1+exp(r(mean)))

tsset datem
gen r3=r1+(rho^1)*F12.r1+(rho^2)*F24.r1
gen r5=r3+(rho^3)*F36.r1+(rho^4)*F48.r1
/* ********************************************************************* */
* Run VAR on RHS in predictive regression, i.e. STGE2
keep if month(dofm(datem))==12
/* ********************************************************************* */
tsset datem
quietly: regress STGE2 L12.STGE2
predict v if e(sample), residual

matrix eqn1=e(b)
scalar p=eqn1[1,1]
display p
matrix drop eqn1
* *****************
regress r1 L12.STGE2
predict u if e(sample), residual

matrix eqn2=e(b)
scalar b=eqn2[1,1]
display b
matrix drop eqn2
* *****************
pwcorr u v
* *****************
regress u v 
matrix eqn3=e(b)
scalar b_u_v=eqn3[1,1]
display b_u_v
drop u v
matrix drop eqn3
* *****************
scalar pA=(_N*p+1)/(_N-3)
display "pA=" pA " p="p

scalar bA = b + b_u_v * (1+3*pA)/_N
display "bA=" bA " b=" b

quietly: su L12.STGE2 if r1<.
scalar  LSTGE2_mean=r(mean)
display LSTGE2_mean

quietly: su STGE2 if r1<.
scalar STGE2_mean=r(mean)

quietly: su r1
scalar ret_mean=r(mean)

scalar a = ret_mean - bA * LSTGE2_mean
scalar c = STGE2_mean - pA * LSTGE2_mean

gen u_true= r1  - (a + bA * L12.STGE2)
gen v_true= STGE2 - (c + pA * L12.STGE2)

* sanity check:  errors must average 0
keep if u_true+v_true<.
su *_true
* *****************
count
* *****************
quietly: regress r1  L12.STGE2
matrix eqn1=e(b)
scalar b1=eqn1[1,1]
matrix drop eqn1

quietly: regress r3 L12.STGE2
matrix eqn3=e(b)
scalar b3=eqn3[1,1]
matrix drop eqn3

quietly: regress r5 L12.STGE2
matrix eqn5=e(b)
scalar b5=eqn5[1,1]
matrix drop eqn5

display b "  " b1 "  " b3 "  " b5
tsset datem
eststo clear
eststo: quietly: ivreg2 r1 L12.STGE2 if r5<., bw(13) robust small
eststo: quietly: ivreg2 r3 L12.STGE2 if r5<., bw(37) small robust
eststo: quietly: ivreg2 r5 L12.STGE2 if r5<., bw(61) small robust
esttab *, se b(4) scalar(r2_a)
esttab using $Tables/AppendixE1.csv,unstack drop(_cons) nonotes compress b(4) se(4) star(c 0.10 b 0.05 a 0.01) nogaps scalar(r2_a) append

* bH: hypothetical beta
program ks, rclass
	args bH 
	version 13
	gen i = 1 + floor(_N*runiform())

	gen u_sim=u_true[i]	
	gen v_sim=v_true[i]

	generate  xi_sim = STGE2					if _n==1
	replace   xi_sim = c 	+  pA * xi_sim[_n-1] 	+ v_sim		if _n>1
	generate ret_sim=  a	+ `bH'* xi_sim[_n-1]	+ u_sim		

	gen ret_yr1=ret_sim
	gen ret_yr3=ret_yr1+rho^1*ret_yr1[_n+1]+rho^2*ret_yr1[_n+2]
	gen ret_yr5=ret_yr3+rho^3*ret_yr1[_n+3]+rho^4*ret_yr1[_n+4]

	generate Lxi_sim = xi_sim[_n-1] if _n>1
	regress ret_yr1 Lxi_sim if r5<.
	return scalar ols1=_b[Lxi_sim]
	
	regress ret_yr3 Lxi_sim	if r5<.
	return scalar ols3=_b[Lxi_sim]

	regress ret_yr5 Lxi_sim
	return scalar ols5=_b[Lxi_sim]

	drop i *_sim ret_yr* 
end

putexcel set $Tables/AppendixE1_pvalues.xlsx, sheet(STGE2) modify
	putexcel A1=("Slope H0") B1=("Dep Variable") C1=("Slope 1Yr") D1=("p-value 1Yr") E1=("Slope 3Yr") F1=("p-value 3Yr") G1=("Slope 3Yr") H1=("p-value 5Yr")

set seed 123456 
local row=2
local reps=5000
forvalues coef=0(-1)0 {
	preserve
	putexcel A`row'=`coef' 
	putexcel B`row'="STGE2" 
	simulate slope1=r(ols1) slope3=r(ols3) slope5=r(ols5), reps(`reps') nodots: ks `coef' 
* *******
	su slope1,d 
	putexcel C`row'=`r(mean)'
	quietly: count if slope1<b1
	putexcel D`row'=`r(N)'/`reps'
	display `r(N)'/`reps'
* *******
	su slope3,d 
	putexcel E`row'=`r(mean)'
	quietly: count if slope3<b3
	putexcel F`row'=`r(N)'/`reps'
	display `r(N)'/`reps'
* *******
	su slope5,d 
	putexcel G`row'=`r(mean)'
	quietly: count if slope5<b5
	putexcel H`row'=`r(N)'/`reps'
	display `r(N)'/`reps'
* *******
	restore
	local row=`row'+1
}
